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Abstract:  Before  corresponding  points  in  images  taken  with  two  cameras  can  be 
used  to  recover  distances  to  objects  in  a  scene,  one  has  to  determine  the  position  and 
orientation  of  one  camera  relative  to  the  other.  This  is  the  classic  photogrammetric 
problem  of  relative  orientation,  central  to  the  interpretation  of  binocular  stereo 
information.  Iterative  methods  for  determining  relative  orientation  were  developed 
long  ago;  without  them  we  would  not  have  most  of  the  topographic  maps  we  do 
today.  Relative  orientation  is  also  of  importance  in  the  recovery  of  motion  and 
shape  from  an  image  sequence  when  successive  frames  are  widely  separated  in  time. 
Workers  in  motion  vision  are  rediscovering  some  of  the  methods  of  photogram metry. 

Described  here  is  a  particularly  simple  iterative  scheme  for  recovering  relative 
orientation  that,  unlike  existing  methods,  does  not  require  a  good  initial  guess  for 
the  baseline  and  the  rotation.  The  data  required  is  a  set  of  pairs  of  corresponding 
rays  from  the  two  projection  centers  to  points  in  the  scene.  It  is  well  known  that  at 
least  five  pairs  of  rays  are  needed.  Less  appears  to  be  known  about  the  existence  of 
multiple  solutions  and  their  interpretation.  These  issues  are  discussed  here  in  detail. 
The  unambiguous  determination  of  all  of  the  parameters  of  relative  orientation  is 
not  possible  when  the  observed  points  lie  on  a  critical  surface.  These  surfaces  and 
their  degenerate  forms  are  analysed  here  as  well.  — 
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to  recover  distances  to  objects  in  a  scene,  one  has  to  determine  the  position 
and  orientation  of  one  camera  relative  to  the  other.  This  is  the  classic 
photograrametric  problem  of  "relative  orientation,"  central  to  the  interpretatic n 
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in  the  recovery  of  motion  and  shape  from  an  image  sequence  when  successive 
frames  are  widely  separated  in  time.  Workers  in  motion  vision  are  redis¬ 
covering  some  of  the  methods  of  photogrammetry. 

Described  here  is  a  particularly  simple  iterative  scheme  for  recovering 
relative  orientation  that,  unlike  existing  methods,  does  not  require  a  good 
initial  guess  for  the  baseline  and  the  rotation.  The  data  required  is  a  set 
of  pairs  of  corresponding  rays  from  the  two  projection  centers  to  points 
in  the  scene.  It  is  well  known  that  at  least  five  pairs  of  rays  are  needed. 
Less  appears  to  be  known  about  the  existence  of  multiple  solutions  and  their 
interpretation.  These  issues  are  discussed  here  in  detail.  The  unambiguous 
determination  of  all  of  the  parameters  of  relative  orientation  is  not 
possible  when  the  observed  points  lie  on  a  critical  surface.  These  surfaces 
and  their  degenerate  forms  are  analysed  here  as  well. 
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1.  Introduction 

The  positions  of  corresponding  points  in  two  images  can  be  used  to  determine 
the  positions  of  points  in  the  environment,  provided  that  the  position  and 
orientation  of  one  camera  with  respect  to  the  other  is  known.  Given  the 
internal  geometry  of  the  camera,  including  its  focal  length  and  the  location  of 
the  principal  point,  rays  can  be  constructed  by  connecting  the  points  in  the 
images  to  their  corresponding  projection  centers.  These  rays,  when  extended, 
intersect  at  the  point  in  the  scene  that  gave  rise  to  the  image  points.  This  is 
how  binocular  stereo  data  is  used  to  determine  the  positions  of  point  s  in  the 
environment  after  the  correspondence  problem  has  been  solved. 

It  is  also  the  method  used  in  motion  vision  when  feature  points  are  tracked 
and  the  image  displacements  that  occur  in  the  time  between  two  successive 
frames  are  relatively  large  (see  for  example  [Ullman  1979]  and  [Tsai  &:  Huang 
1984]).  The  connection  between  these  two  problems  has  not  attracted  much 
attention  before,  nor  has  the  relationship  of  motion  vision  to  some  aspects 
of  photogrammetry  (but  see  [Longuet- Higgins  1981]).  It  turns  out,  for  exam¬ 
ple,  that  the  well  known  motion  field  equations  [Longuet-Higgins  &  Prazdny 
1980.  Bruss  &  Horn  1983]  are  just  the  parallax  equations  of  photogrammetry 
[Hallert  1960,  Moffit  &  Mikhail  1980]  that  occur  in  the  incremental  adjust¬ 
ment  of  relative  orientation.  Most  papers  on  relative  orientation  only  give  the 
equation  for  y-parallax,  corresponding  to  the  equation  for  the  y-component  of 
the  motion  field  (see  for  example  the  first  equation  in  [Gill  1964],  equation  (1) 
in  [Jochmann  1965],  and  equation  (6)  in  [Oswal  1967]).  Some  papers  actually 
give  equations  for  both  x-  and  y-parallax  (see  for  example  equation  (9)  in 
[Bender  1967]). 

In  both  binocular  stereo  and  large  displacement  motion  vision  analysis, 
it  is  necessary  to  first  determine  the  relative  orientation  of  one  camera  with 
respect  to  the  other.  The  relative  orientation  can  be  found  if  a  sufficiently  large 
set  of  [)airs  of  corresponding  image  points  have  been  identified  [Thompson 
1959b,  Thompson  1968,  Ghosh  1972,  Schwidefsky  1973,  Slama  et  al.  1980, 
Moffit  k  Mikhail  1980,  Wolf  1983,  Horn  1986]. 

Let  us  use  the  terms  left  and  right  to  distinguish  the  two  cameras  (in  the 
ca.se  of  the  application  to  motion  vision  these  will  be  the  camera  positions  and 
orientations  corresponding  to  the  earlier  and  the  later  frames  respectively).^ 
The  ray  from  the  center  of  projection  of  the  left  camera  to  the  center  of 
projection  of  the  right  camera  is  called  the  baseline  (see  Fig.  1).  A  coordinate 
system  can  be  erected  at  each  projection  center,  with  one  axis  along  the 

’ll!  what  follows  we  use  the  coorilinale  system  of  the  right  (or  later)  camera  as 
the  reference.  One  can  simply  interchange  left  and  right  if  it  happens  to  be  more 
convenient  to  use  the  coordinate  system  of  the  left  (or  earlier)  camera. 
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optical  axis,  that  is,  perpencliciil2u:  to  the  image  plane.  The  other  two  axes 
are  in  each  case  parallel  to  two  convenient  orthogonal  directions  in  the  image 
plane  (such  as  the  edges  of  the  image,  or  lines  connecting  pairs  of  fiducial 
marks).  The  rotation  of  the  left  camera  coordinate  system  with  respect  to 
the  right  is  called  the  orientation. 
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Figure  1.  Points  in  the  environment  are  viewed  from  two  camera  positions. 
The  relative  orientation  is  the  direction  of  the  baseline  b,  and  the  rotation 
R  relating  the  left  and  right  coordinate  systems.  The  directions  of  rays 
to  at  least  five  scene  points  must  be  known  in  both  camera  coordinate 
systems. 


Note  that  we  cannot  determine  the  length  of  the  baseline  without  knowledge 
about  the  length  of  a  line  in  the  scene,  since  the  ray  directions  are  unchanged 
if  we  scale  all  of  the  distances  in  the  scene  and  the  baseline  by  the  same 
positive  sczde-factor.  This  means  that  we  should  treat  the  baseline  as  a  unit 
vector,  and  that  there  are  really  only  five  unknowns — -three  for  the  rotation 
^lnd  two  for  the  direction  of  the  baseline.^ 

2.  Existing  Solution  Methods 

Various  empirical  procedures  have  been  devised  for  determining  the  relative 
orientation  in  an  analog  fashion.  Most  commonly  used  are  stereoplotters,  op- 

*If  we  treat  the  baseline  as  a  unit  vector,  its  actual  length  becomes  the  unit  of 
length  for  all  other  quantities. 
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tical  devices  that  permit  viewing  of  image  pairs  and  superimposed  synthetic 
features  called  floating  marks.  Differences  in  ray  direction  parallel  to  the 
baseline  are  called  horizontal  disparities  (or  i-parallaxes),  while  differences  in 
ray  direction  orthogonal  to  the  baseline  are  called  vertical  disparities  (or  y- 
parcdlaxes).*  Horizontal  disparities  encode  distances  to  points  on  the  surface 
and  axe  the  quantities  sought  after  in  measurement  of  the  underlying  topogra¬ 
phy.  There  .should  be  no  vertical  disparities  when  the  device  is  adjusted  to  the 
correct  relative  orientation,  since  the  rays  from  the  left  and  right  projection 
center  must  lie  in  a  plane  that  contains  the  baseline  if  they  axe  to  intersect. 

The  methods  used  in  practice  to  determine  the  correct  relative  orien¬ 
tation  depend  on  successive  adjustments  to  eliminate  the  vertical  disparity 
at  each  of  five  or  six  carefully  chosen  points  [Sailor  1960,  Thompson  1964, 
Slama  et  al.  1980,  Moffit  Sz  Mikhail  1980,  Wolf  1974].  In  each  of  the  ad¬ 
justments  a  single  parameter  of  the  relative  orientation  is  varied  in  order  to 
remove  the  vertical  disparity  at  one  of  the  points.  Which  adjustment  is  made 
to  eliminate  the  vertical  disparity  at  a  particular  point  depends  on  the;  partic¬ 
ular  method  is  chosen.  In  each  case,  however,  one  of  the  adjustments,  rather 
than  being  guided  visually,  is  made  by  an  eimount  that  is  calculated,  using 
the  measured  values  of  earlier  adjustments.  The  calculation  is  based  on  the 
assumptions  that  the  surface  being  viewed  can  be  approximated  by  a  plane, 
that  the  baiseline  is  roughly  parallel  to  this  plane,  and  that  the  optical  axes  of 
the  two  cameras  axe  roughly  perpendicular  to  this  plane.  The  whole  process 
is  iterative  in  nature,  since  the  reduction  of  vertical  disparity  at  one  point 
by  means  of  an  adjustment  of  a  single  parameter  of  the  relative  orientation 
disturbs  the  vertical  disparity  at  the  other  points.  Convergence  is  usually 
rapid  if  a  good  initial  guess  is  available.  It  can  be  slow,  however,  when  the 
assumptions  on  which  the  calculation  is  based  axe  violated,  such  as  in  “acci- 
dented”  or  hilly  terrain  [Van  Der  Weele  1959-60).  These  methods  typically 
use  Euler  angles  to  represent  three-dimensional  rotations  [Korn  &  Korn  1968] 
(traditionally  denoted  by  the  greek  letters  k,  4>,  and  w).  Euler  angles  have 
a  number  of  short-comings  for  describing  rotations  that  become  particularly 
noticable  when  these  angles  become  large. 

There  also  exist  related  digital  procedures  that  converge  rapidly  when  a 
good  initial  guess  of  the  relative  orientation  is  available,  as  is  usually  the  case 
when  one  is  interpreting  aerial  photography  [Slama  et  al.  1980].  Not  all  of 
these  methods  use  Euler  angles.  Thompson  [1959b],  for  example,  uses  twice 
the  Gibb’s  vector  [Korn  &  Korn  1968]  to  represent  rotations.  These  proce- 

^This  naming  convention  stems  from  the  observation  that,  roughly  speaking,  in 
the  usual  viewing  arrangement,  horizontal  disparities  correspond  to  left- right 
di.splaremcnts  in  the  image,  whereas  vertical  disparities  correspond  to  up-down 
displacements. 


dures  may  fail  to  converge  to  the  correct  sohition  when  the  initial  guess  is  far 
off  the  mark.  In  the  application  to  motion  vision,  approximate  translational 
and  rotational  components  of  the  motion  are  often  not  known  initially,  so  a 
procedure  that  depends  on  good  initial  guesses  is  not  particulaily  useful. 

Normally,  the  directions  of  the  rays  are  obtained  from  points  in  images 
generated  by  projection  onto  a  planar  surface.  In  this  case  the  directions  are 
confined  to  the  field  of  view  as  determined  by  the  active  area  of  the  image  plane 
and  the  distance  to  the  center  of  projection.  The  field  of  view  is  always  less 
than  a  hemi-sphere,  since  only  points  in  front  of  the  camera  can  be  imaged.'’ 
The  method  described  here  applies,  however,  no  matter  how  the  directions 
to  points  in  the  scene  are  determined.  There  is  no  restriction  on  the  possible 
directions  of  the  rays.  We  do  assume,  however,  that  we  can  tell  which  of  two 
opposite  semi-infinite  line-segments  the  point  lies  on.  If  a  point  lies  on  the 
correct  line-segment  we  will  say  that  it  lies  in  front  of  the  camera,  otherwise 
it  will  be  considered  to  be  behind  the  camera. 

The  problem  of  relative  orientation  is  generally  considered  solved,  and 
so  has  received  little  attention  in  the  photogrammetric  literature  in  recent 
times  [Van  Der  Weele  1959-60].  In  the  annual  index  of  Photogrammetric 
Engineering,  for  example,  there  is  only  one  reference  to  the  subject  in  the  last 
ten  years  [Ghilani  1983]  and  six  in  the  decade  before  that.  This  is  very  little 
in  comparison  to  the  large  number  of  papers  on  this  subject  in  the  fifties, 
3is  well  as  the  sixties,  including  [Gill  1964],  [Sailor  1965],  [Jochmeinn  1965], 
[Ghosh  1966],  [Forrest  1966]  and  [Oswal  1967]. 

In  this  paper  we  discuss  the  relationship  of  relative  orientation  to  the 
problem  of  motion  vision  in  the  situation  where  the  motion  between  the  ex¬ 
posure  of  successive  frames  is  relatively  large.  Also,  a  new  iterative  algorithm 
is  described  here,  as  well  as  a  way  of  dealing  with  the  situation  when  there 
is  no  initial  guess  available  for  the  rotation  or  the  direction  of  the  baseline. 
The  advantages  of  the  unit  quaternion  notation  for  representing  rotations  are 
illustrated  as  well.  FineUy,  we  discuss  critical  surfaces,  surface  shapes  that 
lead  to  difficulties  in  establishing  a  unique  relative  orientation. 

3.  Coplanarity  Condition 

If  the  ray  from  the  left  camera  and  the  corresponding  ray  from  the  right 
camera  are  to  intersect,  they  must  to  lie  in  a  plane  that  also  contains  the 
baseline.  Thus,  if  b  is  the  vector  respresenting  the  baseline,  Tr  the  ray  from 
the  right  projection  cent(T  to  the  point  in  the  scene  and  r/  the  ray  from  the 

■’The  field  of  view  is  larger  than  a  hemi-sphere  in  some  fish-eye  lenses,  where  there 
is  significant  radial  distortion. 
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left  projection  center  to  the  point  in  the  scene,  then  the  triple  product 

(b  Tr] 

equals  zero,  where  r{  =  Rot(r/)  is  the  left  ray  rotated  into  the  right  camera’s 
coordinate  system.®  This  is  the  coplanarity  condition  (see  Fig.  2). 


\  / 
\  / 


Figure  2.  Two  rays  approach  closest  where  they  are  intersected  by  a 
line  perpendicular  to  both.  If  there  is  no  measurement  error,  and  the 
relative  orientation  has  been  recovered  correctly,  then  the  two  rays  actually 
intersect.  In  this  case  the  two  rays  and  the  baseline  lie  in  a  common  plane. 


We  obtain  one  such  constraint  from  each  pair  of  rays.  There  will  be  an  infinite 
number  of  solutions  for  the  baseline  and  the  rotation  when  there  are  fewer 
than  five  pairs  of  rays,  since  there  are  five  unknowns  and  each  pair  of  rays 
yields  only  one  constraint.  Conversely,  if  there  are  more  than  five  pairs  of 
rays,  the  constraints  are  likely  to  be  inconsistent  as  the  result  of  small  errors 
in  the  measurements.  In  this  case,  no  exact  solution  of  the  set  of  constraint 
equations  exists,  and  it  makes  sense  instead  to  minimize  the  sum  of  squares  of 
(errors  in  the  constraint  equations.  In  practice,  one  should  use  more  than  five 
pairs  of  rays  in  order  to  reduce  the  influence  of  mesisurement  errors  [Jochmann 
1965.  Ghosh  1966].  We  shall  see  later  that  the  added  information  also  allows 
one  to  eliminate  spurious  sohitions. 

®The  baseline  vector  b  is  here  also  measured  in  the  coordinate  system  of  the  right 
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4.  What  is  the  Appropriate  Error  Term? 

The  triple  product  [b  rj  rr)  is  zero  when  the  left  and  right  ray  are  coplanar 
with  the  baseline.  It  is  not  immediately  apparent,  however,  that  the  triple 
product  itself  is  necess^lrily  the  ideal  measure  of  departure  from  best  fit.  It  is 
worthwhile  exploring  the  geometry  of  the  two  rays  more  carefully.  Consider 
the  points  on  the  rays  where  they  approach  each  other  the  closest  (see  Fig,  2). 
The  line  connecting  these  points  will  be  perpendicular  to  both  rays,  and  hence 
parallel  to  (rj  x  r^).  As  a  consequence,  we  can  write 

ar;  +  7(r;xr^)  =  b  +  ^rp, 

where  a  and  /3  are  proportional  to  the  distances  along  the  left  and  the  right  ray 
to  the  points  where  they  approach  most  closely,  while  7  is  proportional  to  the 
shortest  distance  between  the  rays.  We  can  find  7  by  taking  the  dot-product 
of  the  equality  above  with  r\  x  Tr.  We  obtain 

7  ||r;  X  Trf  =  [b  r)  r^]. 

Similarly,  taking  the  dot-products  with  Fr  x  (r{  x  Fr)  and  fJ  x  (f}  x  Fr),  we 
obtain 

O  \\r'i  X  Frll^  =  (b  X  Fr)  •  (f',  X  Fr), 

/?||f;  X  Frf  =  (b  X  f',)  •  (f)  X  Fr). 

The  magnitudes  of  o  and  /?  are  the  distances  along  the  rays  to  the  point  of 
closest  approach  when  Fr  and  fJ  are  unit  vectors.  It  turns  out,  however,  that 
we  are  more  interested  here  in  the  signs  than  in  the  magnitudes  of  a  and  /3. 

Normally,  the  points  where  the  two  rays  approach  the  closest  will  be  in 
front  of  both  cameras,  that  is,  both  a  and  will  be  positive.  If  the  estimated 
bziseline  or  rotation  is  in  error,  however,  then  it  is  possible  for  one  or  both  of 
the  calculated  parameters  a  and  /?  to  be  negative.  We  will  use  this  observation 
later  to  distinguish  amongst  different  apparent  solutions. 

We  have  shown  that  the  perpendicular  distance  between  the  left  and  the 
right  ray  is  equal  to  ratio  of  the  triple  product  [b  r}  Fr]  to  the  magnitude 
squared  of  (r|  x  Fr).  But  the  measurement  errors  are  in  the  image,  not  in 
the  scene.  Thus  a  least-squares  procedure  should  be  based  on  an  error  in 
determining  the  direction  of  the  rays,  not  on  the  distance  of  closest  approach. 
To  arrive  at  such  a  measure,  we  can  look  at  the  angle,  9  say,  between  the 
projections  of  the  left  and  right  ray  into  a  plane  perpendicular  to  b  (see 
Fig.  3).  This  angle  will  be  zero  when  the  vertical  disparity  is  zero,  that  is, 
v’  en  the  left  and  right  rays  are  coplanar  with  the  baseline. 

The  projections  of  r)  and  Fr  into  a  plane  perpendicular  to  b  are  given  by 

r'f  =  f'i  —  (f)  •  b)b  =  (b  X  f'/)  X  b, 
fr  =  Fr  -  (Fr  b)b  =  (b  X  Fr)  X  b, 


J 


Relative  Orientation 


7 
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Figure  3.  One  measure  of  the  departure  from  the  copianarity  condition 
is  the  angle  6  between  the  planes  formed  by  the  left  ray  and  the  baseline 
and  the  plane  formed  by  the  right  ray  and  the  baseline.  The  angle  may 
be  found  by  projection  the  two  rays  onto  a  plane  perpendicular  to  the 
baseline. 


where  we  have  used  the  fact  that  b  is  a  unit  vector.  It  is  easy  to  show  that 


||r;||  =  ||b  X  r}||  and  ||r.||  =  ||b  x  r.|| , 


keeping  in  mind  again  that  b  is  a  unit  vector,  and  that  b  x  r}  and  b  x  Tr  axe 
perpendicular  to  b.  The  cross-product  of  the  two  projected  vectors  r|  and  Fr 
will  be  parallel  to  the  baseline  and  have  magnitude  proportional  to  the  sine 
of  the  angle  between  the  projected  vectors,  and  so, 

[b  f;  F^]  [b  r}  Vr] 


sin  B  = 


iFillll^rll  l|b  X  r;i|  l|b  X  r.l 


where  we  have  used  the  fact  that  [b  rj  Tr]  =  [b  r;  Xr],  something  that  can 
easily  be  verified. 


We  could  use  the  sine  of  the  angle  between  the  projected  vectors  directly 
as  a  measure  of  departure  from  best  fit.  This  is  not  as  good  an  idea  as  it  may 
appear  at  first  sight,  because  of  what  happens  when  one  of  the  rays  becomes 
nearly  parallel  to  the  baseline.  In  this  case  the  angle  will  vary  rapidly  with 
small  changes  in  the  direction  of  the  ray.  Correspondingly,  one  of  the  terms 
in  the  denominator  in  the  expression  for  the  sine  of  the  angle  becomes  small. 
It  is  better  to  normalize  the  expression  by  multiplying  by  the  lengths  of  the 
projected  vectors.  Then  we  obtain  the  area  of  the  parallelogram  formed  by 
the  projections  of  the  rays  into  a  plane  perpendicular  to  the  baseline,  namely, 

11*1(1  l|rr(|  sin^  =  [b  f'l  Fr]  =  [b  r{  r^]. 

This  discussion  confirms  that  the  triple  product  itself  is  a  good  measure  of  the 
departure  from  best  fit.  This  is  convenient,  since  it  makes  the  least  squares 
analysis  reasonably  straightforward.  If  one  desires  to  use  a  different  error 
measure,  one  can  weight  the  terms  in  the  sums  to  follow  accordingly. 

5.  Least  Squares  Solution  for  the  Baseline 

If  the  rotation  is  known,  it  is  easy  to  find  the  best  fit  baseline,  as  we  show  next. 
This  is  useful,  despite  the  fact  that  we  do  not  usually  know  the  rotation.  The 
reason  is  that  the  ability  to  find  the  best  baseline,  given  a  rotation,  reduces 
the  dimensionality  of  the  search  space  from  five  to  three. 

Let  {r/_,}  and  {rr,j},  for  i  =  l...n,  be  corresponding  sets  of  left  and 
right  rays.  We  wish  to  minimize 

E  =  ^[b  r;_,  r^,.]^  =  •  (r',,.  x 

1=1  1=1 

subject  to  the  condition  ||b||^  =  1,  where  rj  ■  is  the  rotated  left  ray  r/^,,  as 
before.  If  we  let  c,  =  rj  ^  x  Tr.i,  we  can  rewrite  the  sum  in  the  simpler  form 


^  =  E(bc.)"  =  b^l^c,cr|b, 

i=l  \i=l  / 

where  we  have  used  the  equivalence  b  •  =  b^c^,  which  depends  on  the 

interpretation  of  column  vectors  as  3  x  1  matrices.  The  term  c,c^  is  a  dyadic 
product,  a  3  X  3  matrix  obtained  by  multiplying  a  3  x  1  matrix  by  a  1  x  3 


matrix. 


The  error  sum  is  a  quadratic  form  involving  the  real  symmetric  matrix 

c  =  f^c.cr. 

1=1 

The  minimum  of  such  a  quadratic  form  is  the  smallest  eigenvalue  of  the  matrix 
C,  attained  when  b  is  the  corresponding  unit  eigenvector  (see,  for  example. 
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the  discussion  of  Rayleigh’s  quotient  in  [Korn  Korn  1968]).  This  can  be 
verified  by  introducing  a  Lagraugian  multiplier  A  and  minimizing 

E'  =  b^Cb  +  A(l-b^’b), 

subject  to  the  condition  b^b  =  1.  Differentiating  with  respect  to  b  and 
setting  the  result  equal  to  zero  yields; 

Cb  =  Ab. 

The  error  corresponding  to  a  particular  solution  of  this  equation  is  found  by 
premultiplying  by  b^: 

E  =  h^'C  b  =  Ab^b  =  A. 

The  three  eigenvalues  of  the  real  symmetric  matrix  C  are  non-negative,  anrl 
can  be  found  in  closed  form  by  solving  a  cubic  equation,  while  each  of  the 
corresponding  eigenvectors  has  components  that  are  the  solution  of  dirce  ho¬ 
mogeneous  equations  in  three  unknowms  [Korn  &:  Korn  196S].  If  the  data 
are  relatively  free  of  measurement  error,  then  the  smallest  eigenvalue  will 
be  much  smaller  than  the  other  two,  and  a  reasonable  approximation  to  the 
sought-after  result  can  be  obtained  by  solving  for  the  eigenvector  using  the 
assumption  that  the  smallest  eigenvalue  is  actually  zero.  This  way  one  need 
not  even  solve  the  cubic  equation  (see  also  [Horn  iz  Weldon  1988]). 

If  b  is  a  unit  eigenvector,  so  is  — b.  Changing  the  sense  of  the  baseline 
does  not  change  the  magnitude  of  the  error  term  [b  rj  Tr].  It  does,  however, 
change  the  signs  of  a,  (3  and  7.  One  can  decide  which  sense  of  the  baseline 
direction  is  appropriate  by  determining  the  signs  of  Oj  and  /?;  for  i  =  1 . . .  n. 
Ideally,  they  should  all  be  positive.  The  solution  for  the  optimal  baseline  is 
not  unique  unless  there  are  at  least  two  pairs  of  corresponding  rays.  The 
reason  is  that  the  eigenvector  we  are  looking  for  is  not  uniquely  determined 
if  more  than  one  of  the  eigenvalues  is  zero,  and  the  matrix  has  rank  less  than 
two  if  it  is  the  sum  of  fewer  than  two  dyadic  products  of  independent  vectors. 
This  is  not  a  significant  restriction,  however,  since  we  need  at  least  five  pairs 
of  rays  to  solve  for  the  rotation  anyway. 

6.  Iterative  Improvement  of  Relative  Orientation. 

If  one  ignores  the  orthonormality  of  the  rotation  matrix,  a  set  of  nine  homoge¬ 
neous  linear  equations  can  be  obtained  by  a  transformation  of  the  coplanarity 
conditions  that  was  first  described  in  [Thompson  1959b].  These  equations 
can  be  solved  when  eight  pairs  of  corresponding  ray  directions  an*  known 
[Longuet-Higgins  1981].  This  is  not  a  least-squares  method  that  can  make 
use  of  redundant  measurements,  nor  can  it  be  applied  when  fewer  than  eight 
points  are  given.  The  method  is  also  strongly  affected  by  measurement  errors 
and  fails  ff)r  certain  configurations  of  points  [Longuet-Higgins  1984]. 
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No  true  closed-form  solution  of  the  least-squares  problem  has  been  found 
for  the  general  case,  where  both  baseline  and  rotation  are  unknown.  However, 
it  is  possible  to  determine  how  the  overall  error  is  affected  by  small  changes 
in  ^he  baseline  and  small  changes  in  the  rotation.  This  allows  one  to  make 
iterative  adjustments  to  the  baseline  and  the  rotation  that  reduce  the  sum  of 
squares  of  errors. 

We  czm  represent  a  small  change  in  the  baseline  by  an  infinitesimal  quan¬ 
tity  6b.  If  this  change  is  to  leave  the  length  of  the  baseline  unaltered,  then 

\\h  +  6hf  =  \\h\\\ 

or 

||bf -h2b.6b  +  ||6bf  =  ||b|^ 

If  we  ignore  quantities  of  second-order,  we  obtain 

6b  ■  b  =  0, 

that  is,  6b  must  be  perpendicular  to  b. 

A  small  change  in  the  rotation  can  be  represented  by  an  infinitesimal 
rotation  vector  Su.  The  direction  of  this  vector  is  parallel  to  the  axis  of  rota¬ 
tion,  while  its  magnitude  is  the  angle  of  rotation.  This  incremental  rotation 
takes  the  rotated  left  ray,  r{,  into 

r'l  =  r}  +  (6uf  X  r}). 

This  follows  from  Rodrigues'  formula  for  the  rotation  of  a  vector  r: 

COS0  r  -f-  sin^(w  x  r)  -f  (1  —  cos0)(u>  •  r)i*;, 

if  we  let  =  ||6u;||,  u;  =  6u;/ l|6ti;||,  and  note  that  6a;  is  an  infinitesimal 
quantity.  Finally  then,  we  see  that  [b  rj  Tr]  becomes 

[(b  -I-  6b)  (r{  -f-  6a;  X  r})  Tr], 


or, 

[b  r|  Fr]  +  [6b  r{  Fr]  +  [b  (6a;  x  f})  Fr], 
if  we  ignore  the  term  [6b  (6a;  x  rj)  Fr],  because  it  contains  the  product  of 
two  infinitesimzd  quantities.  We  can  expand  two  of  the  triple  products  in  the 
expression  above  and  obtain 

[b  f',  Fr]  -I-  (f',  X  Fr)  •  6b  +  (Fr  X  b)  •  (6a;  X  f|). 


or 


<  -t-  c  •  6b  -f  d  •  6a;, 


for  short,  where 

<  =  [b  f',  Fr],  c  =  f}  X  Fr,  and  d  =  Fj  X  (Fr  X  b). 
Now,  we  are  trying  to  minimize 


E  =  ^  (t,  -t-  c,  •  6b  -f  d,  •  6a;)^  , 

i=l 
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subject  to  the  condition  b  •  =  0.  This  constrained  minimization  problem 

can  be  transformed  into  an  equivalent  unconstrained  form  by  introduction 
of  a  Lagrangian  multiplier.  Instead  of  minimizing  E  itself,  we  then  have  to 
minimize: 

E' =  £  +  A(b  •  .5b). 

Differentiating  E'  with  respect  to  ^b,  and  setting  the  result  equal  to  zero 
yields, 

n 

(ti  +  c,  •  6b  +  d,  •  6u;)  c,  +  Ab  =  0. 

«=i 

Before  we  can  proceed,  we  need  to  eliminate  the  unknown  Lagrangian  multi¬ 
plier  A.  The  dot-product  of  the  expression  with  b  leads  to 

n 

(ti  -(-  c,  •  6b  -f  d,  •  Su)  c,  •  b  -F  A(b  •  b)  =  0, 

1=1 

which,  since  b  •  b  =  1,  gives  us  a  value  for  A  that  we  can  use  to  compute  the 
term 

n 

Ab  =  -b  (ti  +  Ci  ■  Sb  +  di  ■  6ai)  c,  •  b, 

i=l 

or 

n 

Ab  =  -(bb^)^(fi  +  c<-6b-f  d.  -6u;)c,. 

1=1 

Finally,  substituting  for  Ab  in  the  equation  above  we  obtain 

n 

B  (ti  -t-  c,  •  6b  -f  di  •  Suf)  c,  =  0, 

1=1 

where 

B  =  (I-  bb^) 

is  the  projection  operator  that  removes  components  of  vectors  parallel  to  the 
baseline,  and  I  is  just  the  3x3  identity  matrix.  We  can  conclude  that  the 
equation  above  relates  quantities  in  a  plane  perpendicular  to  the  baseline. 

Finally,  if  we  differentiate  E'  with  respect  to  6a;  and  set  this  result  also 
equal  to  zero,  we  obtain 

n 

(fi  -1-  Ci  •  6b  -I-  d,  ■  6a;)  d,  =0. 

i=l 

Together,  the  two  vector  equations  constitute  six  linear  scalar  equations  in 
the  six  unknown  components  of  6b  and  6a;.  We  can  rewrite  them  in  the  more 
compact  form: 


BC  6b  +  BE  6u>  =  —Be 
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or 


where 


BC  BF\  (  8h\  _  Bc\ 
F'^  D  jysu,)  ~  \  d  ) 


C  =  ^c,c'[,  F^^adf, 


i=l 


i=l 


and  D  =  ^  d,df , 

i=l 


while 

n  n 

c  =  ^  tiCi  and  d  =  ^  <,d,. 
i=l  1=1 

The  matrix  B  is  singular,  as  can  be  seen  by  noting  that  b  is  an  eigenvector  with 
zero  eigenvalue.  Thus  the  first  three  equations  above  are  not  independent. 
One  of  them  will  have  to  be  removed.  Fortunately,  we  also  still  have  to 
incorporate  the  condition  that  be  perpendicular  to  b.  We  czm  do  this  by 
replacing  one  of  the  first  three  equations  with  the  linear  equation  b  ■  6b  =  0. 
For  the  best  numerical  accuracy  one  should  eliminate  the  equation  with  the 
smallest  coefficients. 


The  above  gives  us  a  way  of  finding  small  changes  in  the  baseline  and 
rotation  that  reduce  the  overall  error  sum.  This  method  can  be  applied  itera¬ 
tively  to  locate  a  minimum.  Numerical  experiments  confirm  that  it  converges 
rapidly  when  a  good  initial  guess  is  available.  Incremental  adjustments  can¬ 
not  be  computed  if  the  six-by-six  coefficient  matrix  becomes  singulair.  This 
will  happen  when  there  are  fewer  than  five  pairs  of  rays,  and  for  certain  rare 
configurations  of  points  in  the  scene  (see  the  discussion  of  critical  surfaces 
later  on). 


7.  Adjusting  the  Baseline  and  the  Rotation 

The  iterative  adjustment  of  the  baseline  is  straightforward: 

b"+*  =  b"  +  6b", 

where  b"  is  the  baseline  estimate  at  the  beginning  of  the  n-th  iteration,  while 
6b”  is  the  adjustment  computed  during  the  n-th  iteration,  as  discussed  in 
the  previous  section.  If  6b”  is  not  infinitesimal,  the  result  will  not  be  a  unit 
vector.  We  can,  and  should,  normalize  the  result  by  dividing  by  its  magnitude. 


7.1.  Adjustment  of  Rotation  using  Unit  Quaternions 

Adjusting  the  rotation  is  a  little  harder.  Rotations  are  conveniently  repre¬ 
sented  by  unit  quaternions  [Stuelpnagle  1964,  Salamin  1979,  Taylor  1982, 
Horn  1986,  Horn  1987a].  The  groundwork  for  the  application  of  the  unit 
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quaternion  notation  in  photogrammetry  was  laid  by  Thompson  [1959a]  and 
Schut  [1958-59].  A  positive  rotation  about  the  axis  u  through  an  angle  9  is 
represented  by  the  unit  quaternion 

q  =  cos{B/2)  +  sin(^/2)  uj, 

where  ui  is  assumed  to  be  a  unit  vector.  Composition  of  rotations  corresponds 
to  multiplication  of  the  corresponding  unit  quaternions.  The  rotated  version 
of  a  vector  r  is  computed  using 

c  f  o  o  o  « 

r  =  qrq  , 

where  q*  is  the  conjugate  of  the  quaternion  q,  that  is,  the  quaternion  obtained 
by  changing  the  sign  of  the  vector  part.  Here,  f  is  a  purely  imaginary  quater¬ 
nion  with  vector  part  r,  while  f'  is  a  purely  imaginary  quaternion  with  vector 
part  r'. 

The  infinitesimal  rotation  8lo  corresponds  to  the  quaternion 

=  1  -|-  Su). 

We  can  adjust  the  rotation  q  by  premultiplying  with  8u>,  that  is, 

q"+*  =  8uj^  q". 

If  8u3”'  is  not  infinitesimal,  will  not  be  a  unit  quaternion,  and  so  the  result 
of  the  adjustment  will  not  be  a  unit  quaternion  either.  This  undesirable  state 
of  affairs  can  be  avoided  by  using  either  of  the  two  unit  quaternions 

8uj  =  y^l  —  ||^<*'||^  +  8u, 

or  _ 

=  (1  -|-  8u})l ^1  -f  )]<5u;]]^. 

Alternatively,  one  can  simply  normalize  the  product  by  dividing  by  its  mag¬ 
nitude. 


7.2.  Adjustment  of  Rotation  using  Orthonormal  Matrices 


The  adjustment  of  rotation  is  a  little  trickier  if  orthonormal  matrices  are  used 
to  represent  rotations.  We  can  write  the  relationship 

r'  =  r  -I-  {8u)  X  r). 


in  the  form 


r'  =  r  -f  W  r, 


where  the  skew-symmetric  matrix  W  is  defined  by 

/  0  —8u}z  8ufy  \ 


in  terms  of  the  components  of  rotation  vector  6u;  =  (6ux,Su)j,^Su>z)^.  Conse¬ 
quently  we  may  write  r'  =  Qr,  where  Q  =  I  W,  or 

(1  —Suiz  \ 

S(jJz  1  —6u}z  I  , 

—  Sitly  SiVz  1  / 

One  couid  then  attempt  to  adjust  the  rotation  by  multiplication  of  the  ma¬ 
trices  Q  and  R  as  folows: 

The  problem  is  that  Q  is  not  orthonormal  unless  8u)  is  infinitesimal.  In  prac¬ 
tice  this  means  that  the  rotation  matrix  will  depart  more  and  more  from 
orthonormality  as  more  and  more  iterative  adjustments  are  made.  It  is  pos¬ 
sible  to  re-normalize  this  matrix  by  finding  the  nearest  orthonormal  matrix, 
but  this  is  complicated,  involving  the  determination  of  the  square-root  of  a 
symmetric  matrix  [Horn,  Hilden  &  Negahdaripour  1988].  This  is  one  place 
where  the  unit  quaternion  representation  has  a  distinct  advantage:  it  is  trivial 
to  find  the  nearest  unit  quaternion  to  a  quaternion  that  does  not  have  unit 
magnitude. 

To  avoid  this  problem,  we  should  really  start  with  an  orthonormal  ma¬ 
trix  to  represent  the  incremental  rotation.  We  cau  use  either  of  the  unit 
quaternions 

8Cj  =  \J^  —  -|-  8ui^ 

or 

6a;  =  (1  -f  8Lj)lyJ\  -t-  |16u;|l^, 
to  construct  the  corresponding  orthonormal  matrix 

(9o  +  -  9^  2(9i9j,  -  q^qz)  2{qxqz  +  go9»)  \ 

2(9j,9x  +  qoqz)  ql-ql  +  q]-  qi  ^qyqz  -  qoqx)  j , 
2(9x9*  -  9o9j/)  ^(qzqy+qoqx)  9o-9r-9y  +  9x/ 
where  is  the  scalar  part  of  the  quaternion  6a;,  while  qx,  qy,  qz  are  the  com¬ 
ponents  of  the  vector  part.  Then  the  adjustment  of  rotation  is  accomphshed 
using 

Note,  however,  that  the  resulting  matrices  will  still  tend  to  depart  somewhat 
from  orthonormality  due  to  numerical  inaccuracies.  This  may  be  a  problem 
if  many  iterations  are  required. 


8.  Inherent  Ambiguities 

The  iterative  adjustment  described  above  may  arrive  at  a  number  of  appar¬ 
ently  different  solutions.  Some  of  these  are  just  different  representations  of 
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the  same  solution,  while  others  are  related  to  the  correct  solution  by  a  simple 
transformation.  First  of  all,  ne'e  that  —  q  represents  the  same  rotation  as  q, 
since 

(-q)f(-q*)  =  qfq*. 

That  is,  antipodal  points  on  the  unit  sphere  in  four  dimensions  represent  the 
same  rotation.  We  can  prevent  any  confusion  by  ensuring  that  the  first  non¬ 
zero  component  of  the  resulting  unit  quaternion  is  positive,  or  that  the  largest 
component  is  positive. 

Next,  note  that  the  triple  product,  [b  rj  Tr],  changes  sign,  but  not  mag¬ 
nitude,  when  we  replace  b  with  — b.  Thus  the  two  possible  senses  of  the 
baseline  yield  the  same  sum  of  squares  of  errors.  However,  changing  the  sign 
of  b  does  change  the  signs  of  both  a  and  All  scene  points  imaged  are  in 
front  of  the  camera,  so  the  distances  should  all  be  positive.  In  the  jjresence 
of  noise,  it  is  possible  that  some  of  the  distances  turn  out  to  be  negative,  but 
with  reasonable  data  almost  all  of  them  should  be  positive.  This  allows  us  to 
easily  pick  the  correct  sense  for  the  baseline. 

Not  so  obvious  is  another  possibility:  Suppose  we  turn  all  of  the  left  mea¬ 
surements  through  TT  radians  about  the  baseline,  in  addition  to  the  rotation 
already  determined.  That  is,  replace  r)  with 

r'/  =  2(b  •  r;)b  -  r). 

This  follows  from  Rodrigues’  formula  for  the  rotation  of  a  vector  r: 

cos  ^  r  -f  sin  ^  (w  x  r)  -t-  ( 1  —  cos  9)(iv  ■  r)  u;, 
with  $  =  TT  and  u;  =  b.  Then  the  triple  product  [b  rJ  Tr]  turns  into 
2(b  •  rj)[b  b  r^]  -  [b  r}  rr]  =  -[b  r}  Tr]. 

This,  once  again,  reverses  the  sign  of  the  error  term,  but  not  its  magnitude. 
Thus  the  sum  of  squares  of  errors  is  unaltered.  The  signs  of  a  and  ^  are 
affected,  however,  although  this  time  not  in  as  simple  a  way  as  when  the 
sense  of  the  baseline  was  reversed. 

If  [b  r;  Tr]  =  0,  we  find  that  exactly  one  of  cv  and  changes  sign.  This 
can  be  shown  as  follows:  The  triple  product  will  be  zero  when  the  left  and 
right  rays  are  coplanar  with  the  baseline.  In  this  case  we  have  7  =  0,  and  so 

orj  =  b  -1-  /^Fr, 

Taking  the  cross-product  with  b  we  obtain 

o(rJ  X  b)  =  /?(rr  X  b), 

If  we  now  replace  v\  by  r"  =  2(b  •  r))b  —  r|,  we  have  for  the  new  distances  a' 
and  along  the  rays: 

-a'(r',  X  b)  =  /3'(rr  x  b), 

We  conclude  that  the  product  has  sign  opi)osite  to  that  of  the  product 
a/3.  So  if  a  and  /?  are  both  positive,  one  of  a'  or  /?'  must  be  negative. 


In  the  presence  of  measurement  error  the  triple  product  will  not  be  ex¬ 
actly  equal  to  zero.  If  the  rays  are  nearly  coplamar  with  the  baseline,  however, 
we  find  that  one  of  a  and  almost  2ilways  changes  sign.  With  very  poor  data, 
it  is  possible  that  both  change  sign.  (Even  with  totally  random  ray  directions, 
however,  this  only  happens  27.3%  of  the  time,  as  determined  by  Monte  Carlo 
methods).  In  any  case,  we  can  reject  a  solution  in  which  roughly  hzdf  the 
distances  are  negative.  Moreover,  we  can  find  the  correct  solution  directly  by 
introducing  an  additional  rotation  of  tt  radians  about  the  baseline. 

9.  Remaining  Ambiguity 

If  we  take  care  of  the  three  apparent  two-way  ambiguities  discussed  in  the 
previous  section,  we  find  that  in  practice  a  unique  solution  is  found,  provided 
that  a  sufficiently  large  number  of  ray  pairs  are  available.  That  is,  the  method 
converges  to  the  unique  global  minimum  from  every  possibly  starting  point 
in  parameter  space. 

Loceil  minima  in  the  sum  of  squares  of  errors  appear  when  only  a  few 
more  than  the  minimum  of  five  ray  pairs  are  available  (as  is  common  in  prac¬ 
tice).  This  means  that  one  has  to  repeat  the  iteration  with  different  starting 
values  for  the  rotation  in  order  to  locate  the  global  minimum.  A  starting 
value  for  the  baseline  can  be  found  in  each  case  using  the  closed-form  method 
described  in  section  5.  To  search  the  parameter  space  effectively,  one  needs 
a  way  of  efficiently  sampling  the  space  of  rotations.  The  space  of  rotations 
is  isomorphic  to  the  unit  sphere  in  four  dimensions,  with  antipodal  points 
identified.  The  rotation  groups  of  the  regular  polyhedra  provide  convenient 
mezins  of  uniformly  sampling  the  space  of  rotations.  The  group  of  rotations  of 
the  tetrahedron  has  12  elements,  that  of  the  hexahedron  and  the  octahedron 
has  24,  and  that  of  the  icosahedron  and  the  dodecahedron  has  60  (representa¬ 
tions  of  these  groups  are  given  in  Appendix  A  for  convenience).  One  can  use 
these  as  starting  values  for  the  rotation.  Alternatively,  one  can  just  generate 
a  number  of  randomly  placed  points  on  the  unit  sphere  in  four  dimensions  as 
starting  values  for  the  rotation. 

When  there  are  only  five  pairs  of  rays,  the  situation  is  different  again.  In 
this  case,  we  have  five  non-linear  equations  in  five  unknowns  and  so  in  general 
expect  to  find  a  finite  number  of  exact  solutions.  That  is,  it  is  possible  to 
find  baselines  and  rotations  that  satisfy  the  coplanarity  conditions  exactly 
and  reduce  the  sum  of  squares  of  errors  to  zero.  If  we  ignore  the  three  types 
of  ambiguities  discussed  above  (incuding  the  rotation  by  tt  radians  about  the 
baseline),  then  there  arc  generally  four  distinct  sets  of  baselines  and  rotations 
that  yield  an  exact  solution.  Typically  only  one  of  these  yields  positive  signs 
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for  all  of  the  distances.  These  are  empirical  observations;  I  have  not  been  abh' 
to  prove  that  there  are  genen  Uy  four  solutions  that  satisfy  the  coiilanarity 
conditions. 

10.  Summary  of  the  Algorithm 

Consider  first  the  case  where  we  have  an  initial  guess  for  the  rotation.  We 
start  by  finding  the  best-fit  baseline  direction  using  the  closed-form  method 
described  in  section  5.  We  determine  the  correct  sense  of  the  baseline  by 
choosing  the  one  that  makes  most  of  the  signs  of  the  distances  positive.  Then 
we  proceed  as  follows: 

•  For  each  pair  of  corresponding  image  points,  we  compute  rj  ,,  the  left 
ray  direction  r;  ,  rotated  into  the  right  camera  coordinat»*  system, 
using  the  present  guess  for  the  rotation. 

•  We  then  compute  the  cross-product  =  r) ,  x  Tr,,,  the  double  cross- 
product  di  =  r',  ,  X  (Tr.i  X  b)  and  the  triple-product  t,  =  [b  r)  ■  Tr,,]. 

•  We  accumulate  the  dyadic  products  c,c,^,  c,d^  and  d,d/ ,  as  well 
as  the  vectors  and  fid^.  The  totals  of  these  quantities  over  all 
image  point  pairs  give  us  the  matrices  C,  F,  D  and  the  vectors  c 
and  d. 

•  We  can  now  solve  for  the  increment  in  the  baseline  (5b  and  the  in¬ 
crement  in  the  rotation  Soj  using  the  method  derived  in  section  6. 

•  We  adjust  the  baseline  and  the  rotation  using  the  methods  discussed 
in  section  7,  and  recompute  the  sum  of  the  squares  of  the  error  terms. 

The  new  orientation  parameters  are  then  used  in  the  next  iteration  of  the 
above  sequence  of  steps.  As  is  the  case  with  most  iterative  procedures,  it  is 
hard  to  decide  when  to  stop.  Typically,  the  total  error  becomes  small  after  a 
few  iterations  and  no  longer  decreases  at  each  step,  because  of  limited  accuracy 
in  the  arithmetic  operations.  One  can  stop  the  iteration  the  first  time  the  error 
increases.  Alternatively,  one  can  stop  after  either  a  fixed  number  of  iterations 
or  after  the  error  becomes  less  then  some  predetermined  threshold. 

WIk’ii  the  decision  has  been  made  to  stop  the  iteration,  a  check  of  the 
signs  of  the  distances  along  the  rays  is  in  order.  If  most  of  them  are  negative, 
the  baseline  direction  should  be  reversed.  If  neither  sense  of  the  ba-seline 
direction  yields  mostly  positive  distances,  one  needs  to  consider  the  possibility 
of  a  rotation  through  tt  radians  about  the  baseline  b.  If  this  also  yields  mixed 
signs,  one  is  dealing  with  a  local  extremum  of  the  error  function;  something 
that  will  only  happen  if  the  initial  guess  is  in  fact  not  adequate. 

If  an  initial  guess  is  not  availabh*,  one  proceeds  as  follows; 


For  each  rotation  in  the  chosen  group  of  rotations,  perform  the  above 
iteration  to  obtain  a  candidate  baseline  and  rotation. 


•  Choose  the  solution  that  yields  the  smallest  total  error. 


When  there  are  many  pairs  of  rays,  the  iterative  algorithm  will  converge  to 
the  global  minimum  error  solution  from  any  initial  guess  for  the  rotation. 
There  is  no  need  to  sample  the  space  of  rotations  in  this  case.  Also,  instead 
of  sampling  the  space  of  rotations  in  a  systematic  way  using  a  finite  group 
of  rotations,  one  can  generate  points  randomly  distributed  on  the  surface  of 
the  unit  sphere  in  four-dimensional  space.  This  provides  a  simpler  means  of 
generating  initial  guesses,  although  more  initial  guesses  may  have  to  be  tried 
than  when  a  systematic  procedure  is  used  to  sample  the  space  of  rotations. 

The  method  given  minimizes  the  sum  of  the  squares  of  the  triple  products 

[b  r;  Tr]. 


If  desired,  one  can  modify  it  to  use  some  multiple  of  the  triple  product  as  an 
error  term  by  weighting  the  contribution  to  the  overall  sum.  This  can  lead  to  a 
considerably  more  complex  optimization  problem  if  the  weights  depend  on  the 
unknown  baseline  and  the  unknown  rotation.  This  happens,  for  example,  if 
we  try  to  minimize  the  sum  of  squares  of  the  sines  of  the  angles  corresponding 
to  verticjil  disparity: 


sin^  = 


(b  r}  Tr 


l|bxr',||  llbxrrir 

If  we  assume  that  the  weighting  factors  vary  slowly  during  the  iterative  pro¬ 
cess,  however,  we  can  to  use  the  current  estimates  of  the  baseline  and  rotation 
in  computing  the  weighting  factors,  and  not  tcike  into  account  the  small  vari¬ 
ations  in  the  error  sum  due  to  changes  in  the  weighting  factors.  That  is,  when 
taking  derivatives,  the  weighting  factors  are  considered  constant.  This  is  a 
good  approximation  when  the  parameters  vary  slowly,  as  they  will  when  one 
is  close  to  an  extremum. 


11.  Critical  Surfaces 

In  certain  rare  cases,  relative  orientation  cannot  be  recovered  fully,  even  when 
there  are  five  or  more  pairs  of  rays.  Normally,  each  error  term  varies  linearly 
with  distance  in  parameter  space  from  the  location  of  an  extremum,  and  so  the 
sum  of  squares  of  errors  varies  quadratically.  There  are  situations,  however, 
where  the  error  terms  to  not  vary  linearly  with  distance,  but  quadratically  or 
higher  order,  in  certain  special  directions  in  parameter  space.  In  this  case,  the 
sum  of  squares  of  errors  does  not  vary  quadratically  with  distance  from  the 
extremum,  but  as  a  function  of  the  fourth  or  higher  power  of  this  distance. 
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This  makes  it  very  difficult  to  accurately  locate  the  extremum.  In  an  extreme 
situation,  the  total  error  may  riot  vary  at  all  along  some  curve  in  parameter 
space.  In  this  case,  the  total  error  is  unaffected  by  a  small  change  in  the 
rotation,  as  long  as  this  change  is  accompanied  by  a  corresponding  small 
change  in  the  baseline.  There  is  no  localized  extremum  and  consequently 
the  solution  for  relative  orientation  is  not  unique.  It  turns  out  that  this 
problem  arises  only  when  the  observed  scene  points  lie  on  certain  surfaces 
called  Gefdhrliche  Fldchen  or  critical  surfaces  [Brmidenberger  1947,  Hofmann 
1949,  Zeller  1912,  Schwidefsky  1973].  VVe  show  next  that  only  points  on 
certain  hyperboloids  of  one  sheet  and  their  degenerate  forms  can  leatl  to  this 
kind  of  problem  (see  also  [Horn  1987b]). 

We  could  try  to  find  a  direction  of  movement  in  parameter  space  (iSb,  bu.') 
that  leaves  the  total  error  unaffected  when  given  a  particular  surface.  Instead, 
we  will  take  the  critical  direction  of  motion  in  the  parameter  space  as  givi-n, 
and  try  to  find  a  surface  for  which  the  total  error  is  unchanged. 

Let  R  be  a  point  on  the  surface,  measured  in  the  right  camera  coordinate 
system.  Then 

if  Tr  =  R  and  o  r)  =  b  +  R, 

for  some  positive  a  and  j3.  In  the  absence  of  measurement  errors, 

[b  r;r.]=^[b  (b  +  R)  R]  =  0. 

We  noted  earlier  that  when  we  change  the  baseline  and  the  rotation  slightly, 
the  triple  product  [b  r]  Tr]  becomes 

[(b  +  ^b)  (r',  +  bu;  x  r])  Tr], 
or,  if  we  ignore  higher-order  terms, 

[b  r)  Fr]  +  (r]  X  Fr)  •  6b  -I-  (Fr  X  b)  •  (6u;  X  f',). 

The  problem  we  are  focusing  on  here  arises  when  this  error  term  is  unchanged 
for  small  movement  in  some  direction  in  the  parameter  space.  That  is  when 

(f)  X  Fr)  •  6b  -t-  (Fr  X  b)  •  (6u;  X  r'l)  =  0, 

for  some  6b  and  bu.  Introducing  the  coordinates  of  the  imaged  points  we 
obtain; 

^  ^((b  +  R)  X  R)  •  6b  4-  (R  X  b)  •  (6u;  x  (b  -h  R)))  =  0, 
or 

(R  X  b)  •  (6a;  x  R)  -|-  (R  x  b)  •  (6a;  x  b)  -h  [b  R  6b]  =  0. 

If  we  expand  the  first  of  the  dot-products  of  the  cross-products,  we  can  write 
this  equation  in  the  form 

(R  •  b)(6a;  ■  R)  -  (b  •  6a;)(R  •  R)  -f  L  •  R  =  0, 


where 


L  =  £  X  b,  while  £  =  b  x  +  (5b. 

The  expression  on  the  left-hand  side  contains  a  part  that  is  quadratic  in  R 
and  a  part  that  is  linear.  The  expression  is  clearly  quadratic  in  X,  Y,  and 
Z,  the  components  of  the  vector  R  =  {X,Y,  Z)^ .  Thus  a  surface  leading  to 
the  kind  of  problem  described  above  must  be  a  quadric  surface  [Korn  &  Korn 
1968]. 

Note  that  there  is  no  constant  term  in  the  equation  of  the  surface,  so 
R  =  0  satisfies  the  equation.  This  means  that  the  surface  passes  through  the 
right  projection  center.  It  is  easy  to  verify  that  R  =  —  b  satisfies  the  equation 
also,  which  mezms  that  the  surface  passes  through  the  left  projection  center 
as  well.  In  fact,  the  whole  baseline  (and  its  extensions),  R  =  kh,  lies  in  the 
surface.  This  means  that  we  must  be  dealing  with  a  ruled  quadric  surface.  It 
can  consequently  not  be  an  ellipsoid  or  hyperboloid  of  two  sheets,  or  one  of 
their  degenerate  forms.  The  surface  must  be  a  hyperboloid  of  one  sheet,  or 
one  of  its  degenerate  forms.  Additional  information  about  the  properties  of 
these  surfaces  is  given  in  Appendix  B,  while  the  degenerate  forms  are  explored 
in  Appendix  C. 

It  should  be  apparent  that  this  kind  of  ambiguity  is  quite  rare.  This  is 
nevertheless  an  issue  of  practical  importance,  however,  since  the  accuracy  of 
the  solution  is  reduced  if  the  points  lie  near  some  critical  surface.  A  textbook 
case  of  this  occurs  in  aerial  photography  of  a  roughly  U-shaped  valley  taken 
along  a  flight  line  parallel  to  the  axis  of  the  valley  from  a  height  above  the 
valley  floor  approximately  equal  to  the  width  of  the  valley.  In  this  case  the 
baseline  lies  on  a  circular  cylinder  that  also  lies  close  to  the  surface  on  which 
the  imaged  points  lie.  This  means  that  it  is  close  to  one  of  the  degenerate 
forms  of  the  hyperboloid  of  one  sheet  (see  Appendix  C). 

Note  that  hyperboloids  of  one  sheet  and  their  degenerate  forms  are  ex¬ 
actly  the  surfaces  that  lead  to  ambiguity  in  the  case  of  motion  vision.  The 
coordinate  systems  and  symbols  have  been  chosen  here  to  make  the  correspon¬ 
dence  between  the  two  problems  more  obvious.  The  relationship  between  the 
two  situations  is  nevertheless  not  quite  as  transparent  as  I  had  thought  earlier 
[Horn  1987b].  In  the  case  of  the  ambiguity  of  the  motion  field,  for  example,  we 
aire  dealing  with  a  two-way  ambiguity  arising  from  infinitesimal  displacements. 
Here  we  are  dealing  with  an  infinite  number  of  solutions  arising  from  images 
taken  with  cameras  that  have  large  differences  in  position  and  orientation. 
Also  note  that  the  symbol  6u}  stands  for  a  small  change  in  a  finite  rotation 
here,  while  it  refers  to  a  difference  in  instantaneous  rotational  velocities  in  the 
motion  vision  case. 


12.  Conclusions 


Methods  for  recovering  the  relative  t)rientation  of  two  caiiK  ra.s  ar<'  o*'  inipoi 
tance  in  both  binocular  stereo  and  motion  vision.  A  new  it<‘iative  nierhod  for 
finding  the  relative  orientation  has  lieen  descrilx'd  here.  It  can  be  n'<-d  even 
when  there  is  no  initial  guess  available  for  the  rotation  or  the  bas<'!ine.  Tlu' 
new  method  does  not  use  Euler  angles  to  repre.senl  the  orientation. 

When  there  are  many  pairs  of  corresponding  image  points,  the  ter;i;i\(‘ 
method  finds  the  global  minimum  from  any  starting  point  in  i)aiametei  space. 
Local  minima  in  the  sum  of  scjuares  of  errors  occur,  howt’ver.  when  tlieo 
are  relatively  few  pairs  of  corresponding  image  points  availa))le.  Me  hod  for 
efficiently  locating  the  global  minimum  in  this  case  were  discuss(’d.  Wiien  only 
five  pairs  of  corresponding  image  points  are  given,  several  exact  soli  tions  of 
the  coplanarity  equations  can  be  found.  Typically  only  one  of  rhe.o'  yields 
positive  distances  to  the  points  in  the  scene,  however  This  allows  one  to  pick 
the  correct  solution  even  when  there  is  no  initial  guess  available. 

All  of  these  methods  fail  in  the  rare  case  that  the  scene  jroints  lie  on  a 
critical  surface. 
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Appendix  A — Rotation  Groups  of  Regular  Polyhedra 

Each  of  the  rotation  groups  of  the  regular  polyhedra  can  he  generated  from 
two  judiciously  chosen  ('leinents.  For  convenience,  however,  an  explicit  rej)- 
rcsentation  of  all  of  th(‘  elements  of  each  of  tin*  groups  is  given  here.  The 
number  of  diffen-nt  component  values  occuring  in  the  unit  quaternions  rej)- 
re.senting  the  rotations  can  be  kept  low  by  careful  choice  of  the  alignment  (jf 
the  polyhedron  with  the  coordinate  axes.  The  attitudes  of  the  j)olyh('dra  here 
were  selected  to  minimize  the  number  of  different  numcn'ical  vahu's  that  occur 
in  the  components  of  the  unit  quaternions.  .4  different  representation  of  the 
group  is  obtained  if  the  vector  parts  of  each  of  the  unit  quaternions  is  rotated 
in  the  same  way.  This  just  corresponds  to  the  rotation  group  of  thi'  polyhe¬ 
dron  in  a  different  attitude  with  respect  to  the  underlying  coordinate  system. 
This  observation  leads  to  a  convenient  way  of  generating  finer  systematic  sam¬ 
pling  patterns  of  the  space  <T  rotations  than  tlu'  ones  provich'd  directly  by  the 
rotation  group  of  a  regular  polyhedron  in  a  particular  alignment  \vith  the 
coordinate  axes  (see  also  [Brou  1983]). 

The  components  of  the  unit  quaternions  here  may  take  on  the  values  0 
and  1 ,  as  well  as  the  following: 

V^-1  ,1  1  ,  ,  \/5+l 

a  =  - .  6  =  -,  c  =  — and  a  =  - . 

4  2  ^2’  4 

Her('  ar<'  the  unit  (luaternions  for  the  twelve  elements  of  the  rotation 
group  of  the  tetrahedron: 

(1,  0,  0.  0)  (0,  1.  0.  0)  (0.  0,  1.  0)  (0.  0,  0,  1) 

(h,  h,  b.  b)  {b.  b.  b.-b)  (b.  b,-b.  6)  {b,  b,-b.-})) 

{b.-b,  b,  b)  (b.-b.  b.-b)  (b.-b.-b,  b)  (b,-b.-b.-h) 

Here  are  the  unit  quaternions  for  the  twenty-four  eh'inents  of  the  rotation 
group  of  the  octahedron  and  the  h(‘xahedron  (cube): 

(1.  0,  0,  0)  (0.  1,  0,  0)  (0.  0,  1,  0)  (0,  0,  0,  1) 

(0,  (J,  c.  c)  (0,  0,  r, -r)  (0,  c,  0.  c)  (0,  c,  0,-c) 

(0,  c,  c,  0)  (0,  c,-c,  0)  (c,  0,  0,  c)  (c,  0,  0,-c) 

(c,  0,  c,  0)  (c,  0,-c,  0)  (c,  c,  0,  0)  (c,-c,  0,  0) 

(b.  b,  b.  b)  (b.  6,  h,-b)  (6,  6,-6,  6)  (6,  6.-6, -6) 

(b.  b,  b,  b)  (b.-b,  6,-6)  (6, -6, -6,  6)  (6, -6, -6, -6) 


Here  are  the  unit  quaternions  for  the  sixty  elements  of  the  rotation  group  of 
the  icosahedon  and  the  dodectahedron: 
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Remember  that  chzmging  the  signs  of  all  the  components  of  a  unit  quaternion 
does  not  change  the  rotation  that  it  represents. 


Appendix  B — Some  Properties  of  Critical  Surfaces 

In  this  appendix  we  develop  some  more  of  the  properties  of  the  critical  sur¬ 
faces.  The  equation  of  a  critical  surface  can  be  written  in  the  form 

(R  X  b)  •  X  R)  -t-  L  •  R  =  0, 
or 

(R  •  b)(^a;  •  R)  -  (b  •  6ut){K  •  R)  -h  L  •  R  =  0, 

where 

L  =  f  X  b,  while  I  =  h  x  6ut  +  6b. 

It  is  helpful  to  first  establish  some  simple  relationships  between  the  quemtities 
appearing  in  the  formula  above.  We  start  with  the  observations  that  f  •  b  =  0, 
that  t  ■  6u}  —  Sh  •  6u,  and  /  x  6b  =  —(6b  •  6u;)b. 

We  can  also  expand  L  to  yield, 

L  =  6(i;  —  (b  •  6u;)b  -|-  6b  x  b. 

It  follows  that  L  •  b  =  0,  that  L  •  6b  =  6u>  •  6b,  and 

L  X  b  =  -(b  X  6u;  -f  6b)  = 

L  X  6u;  =  {6b  •  6u;)b  -  (b  •  6a;)/. 

We  have  already  established  that  R  =  1;  b  is  an  equation  for  one  of  the  rulings 
passing  through  the  origin.  A  hyperboloid  of  one  sheet  has  two  intersecting 
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families  of  rulings,  so  there  should  be  a  second  ruling  passing  through  tlu- 
origin.  Consider  the  vector  S  defined  by 

S  =  (L  X  Sco)  X  L, 
which  can  be  written  in  the  form 

S  =  (L  •  L)Sc^  -  (L  ■  Scu)L, 
or 

S  =  ((So;  •  6b)£  +  (b  ■  (5a;)(£  •  £)b, 

so  that  S  •  b  =  (b  •  6u})(£  ■  1)  and  S  •  Su  =  (Siv  ■  Sh)^  +  (b  •  8u}Y(l  •  £}. 

If  we  substitute  R  =  A:  S  into  the  formula 

(R  X  b)  •  (Su)  X  R)  +  L  •  R, 
we  obtain  zero,  since  L  •  S  =  0  and 

S  X  b  =  (6u;  •  (5b)L 

is  orthogonal  to 

S  X  6u>  =  — (L  •  8u))'L  X  8ij. 

We  conclude  that  R  =  A:S  is  an  equation  for  the  other  ruling  that  passes 
through  the  right  projection  center. 

There  are  two  families  of  parallel  pleines  that  cut  an  ellipsoid  in  circular 
cross-sections  [Hilbert  &  Cohn-Vossen  1953].  Similarly,  there  are  two  families 
of  parallel  planes  that  cut  a  hyperboloid  of  one  sheet  in  circular  cross-sections. 
One  of  these  families  consists  of  planes  perpendicular  to  the  baseline,  that  is, 
with  common  normal  b.  We  can  see  this  by  substituting  R  •  b  =  fc  in  the 
equation  of  the  critical  surface.  We  obtain 

k{8u}  ■  R)  -  (b  •  6u;)(R  •  R)  +  L  •  R  =  0, 
or 

(b  •  6u;)(R  ■K)-ik8u  +  L)-R  =  0. 

This  is  the  equation  of  a  sphere,  since  the  only  second-order  term  in  R  is  a 
multiple  of 

R-R  =  X^  +  Y^  +  Z\ 

We  can  conclude  that  the  intersection  of  the  critical  surface  and  the  plane  is 
also  the  intersection  of  this  sphere  and  the  plane,  and  so  must  be  a  circle. 
The  same  applies  to  the  intersection  of  the  critical  surface  and  the  family  of 
planes  with  common  normal  8u},  since  we  get 

(b  •  (5a;)(R  •  R)  -  (A:  b  +  L)  •  R  =  0, 

when  we  substitute  R  ■  8u}  =  k  into  the  equation  of  the  critical  surface. 

The  equation  of  the  critical  surface  is  given  in  the  implicit  form  /(R)  = 
0.  The  equation  of  a  tangent  plane  to  such  a  surface  can  be  obtained  by 
differentiating  with  respect  to  R: 

N  =  (R  X  8u})  X  b  -t-  (R  X  b)  X  -f  L 


N  =  (R  •  b)6w  +  (R  •  ^a;)b  -  2(b  •  6u;)R  +  L 
The  tangent  plEine  at  the  origin  has  normal  L.  This  tangent  plane  contains 
the  baseline  (since  L  •  b  =  0),  as  well  as  the  other  ruling  passing  through  +he 
origin  (since  L  •  S  =  0).  Note  that  the  normal  to  the  tangent  plane  is  not 
constant  along  either  of  these  rulings,  as  they  would  be  if  we  were  dealing 
with  a  developable  surface. 

In  the  above  we  have  not  considered  a  large  number  of  degenerate  situ¬ 
ations  that  cam  occur.  The  reader  is  referred  to  Appendix  C  for  a  detailed 
analysis  of  these. 

Appendix  C — Degenerate  Critical  Surfaces 

There  axe  a  number  of  special  alignments  of  the  infinitesimal  change  in  the 
rotation,  Sui,  with  the  baiseline,  b,  and  the  infinitesimal  change  in  the  baseline 
(5b  that  lead  to  degenerate  forms  of  the  hyperboloid  of  one  sheet. 

One  of  the  rulings  passing  through  the  origin  is  given  by  R  =  fc  b,  while 
the  other  is  given  by  R  =  fcS.  If  these  two  rulings  become  parallel,  we  are 
dealing  with  a  degenerate  form  that  has  only  one  set  of  rulings,  that  is  a 
conical  surface.  Now 

S  =  {8uj-  8h)t  +  (b  •  6u){e  •  e)h, 

is  parallel  to  b  only  when  (Su)  •  6b)  =  0,  since  t  is  perpendicular  to  b.  In  this 
case 

6b  •  6a;  =  0  and  6b  •  b  =  0, 

so  6b  =  A:(b  X  6a;)  for  some  constant  k.  Consequently  t  =  {k  +  l)(b  x  6a;), 
The  vertex  of  the  conical  surface  must  lie  on  the  baseline  since  the  baseline 
is  a  ruling,  and  every  ruling  passes  through  the  vertex.  It  can  be  shown  that 
the  vertex  actually  lies  at  R  =  — (fc  +  l)b. 

We  also  know  that  cross-sections  in  planes  perpendicular  to  the  baseliiK' 
are  circles.  This  tells  us  that  we  are  dealing  with  elliptical  cones.  Right 
circular  cones  cannot  be  critical  surfaces.  It  can  be  shown  that  the  main  axis 
of  the  elliptical  cone  lies  in  the  direction  b  -|-  Su. 

A  special  case  of  the  special  case  above  occurs  when 

||b  X  6a;||  =  0, 

that  is  Su  II  b.  Here  6a;  =  fcb  for  some  constant  k  and  so  f  =  6b  and 
L  =  6b  X  b.  The  equation  of  the  surface  becomes 

it(R  •  b)2  -  k{h  ■  b)(R  •  R)  -f  L  ■  R  =  0, 
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This  is  the  equation  of  a  circular  cylinder  with  axis  parallel  to  the  baseline. 
In  essence,  the  vertex  of  the  c(  ue  has  receded  to  infinity  along  the  baseline. 

Another  special  case  arises  when  the  radius  of  the  circular  cross-sections 
with  planes  perpendicular  to  the  baseline  becomes  infinite.  In  this  case  we 
obtain  straight  lines,  and  hence  rulings,  in  these  planes.  The  hyperbolic 
paraboloid  is  the  degenerate  form  that  has  the  property  that  each  of  its  two 
sets  of  rulings  can  be  obtained  by  cutting  the  surface  with  a  set  of  parallc'l 
planes  [Hilbert  &  Cohn-Vossen  1953].  This  happens  when  8u}  is  perpendicular 
to  b,  that  is,  b  •  =  0.  The  equation  of  the  surface  in  this  case  simplifies  to 

(R  ■  b)(6a;  •  R)  +  L  •  R  =  0. 

The  intersection  of  this  surface  with  any  plane  perpendicular  to  the  bascliiu' 
is  a  straight  line.  We  can  show  this  by  substituting 

R  •  b  =  A:, 

into  the  equation  of  the  surface.  We  obtain 

(A:  <^a>  -f  L)  •  R  =  0, 

that  is,  the  equation  of  another  plane.  Now  the  intersection  of  two  planes  is 
a  straight  line.  So  we  may  conclude  that  the  intersection  of  the  surface  and 
the  original  plane  is  a  straight  line.  We  can  show  in  the  same  way  that  the 
intersection  of  the  surface  with  any  plane  perpendicular  to  Su^  is  a  straight 
line  by  substituting 

R  •  =  k 

into  the  equation  of  the  surface.  It  can  be  shown  that  the  saddle  point  of  the 
hyperbolic  paraboloid  surface  lies  on  the  baseline. 

A  special  case  of  particular  interest  arises  when  Su}  is  perpendicular  to 
both  b  and  <5b,  that  is, 

b  •  <!>u;  =  0  and  <5b  •  (5u;  =  0, 

and  so  Su  =  A'((!)b  x  b),  for  some  constant  k.  Then  i  =  (A:  +  l)<5b  and 
L  =  (A:  -f  l)(<^b  X  b).  The  equation  of  the  surface  becomes 

A-(R  •  b)((^b  X  b)  •  R)  -t-  (A-  +  l)(((5b  X  b)  •  R)  =  0, 

or  just 

(A:(R  -  b)  +  (A:  -)-  l))((^b  X  b)  •  R)  =  0. 

so  either 

(6b  X  b)  ■  R  =  0  or  A:(R  •  b)  +  (A:  +  1)  =  0. 

The  first  of  these  is  the  equation  of  a  plane  containing  the  baseline  b  and 
the  vector  6b.  The  second  is  the  equation  of  a  plane  perpendicular  to  the 
baseline.  So  the  solution  degenerates  in  this  case  into  a  surface  consisting  of 
two  int('rs<'cting  planes.  One  of  these  planes  appears  only  as  a  line  in  each  of 
the  two  images,  since  it  passes  through  both  projection  centers,  and  so  does 
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not  really  contribute  to  the  image.  It  is  fortunate  that  planes  can  only  lx- 
degenerate  surfaces  if  they  are  perpendicular  to  the  baseline,  since  surface's 
that  are  almost  planar  occur  frequently  in  aerial  photography.' 

To  summarize  then,  we  have  the  following  degenerate  cases: 


elliptical  cones  when  Su/  ±  Sb, 
circular  cylinders  when  Sof  ||  b, 
hyperbolic  pzu-abloids  when  ±  b, 
intersecting  planes  when  So/  ±  Sb  and  S(ij  ±  b. 


For  further  details,  and  a  proof  that  not  all  hyperboloids  of  one  sheet  passing 
through  the  origin  lead  to  critical  surfaces,  see  [Horn  1987b]. 


'The  baseline  was  nearly  perpendicular  to  the  surface  in  the  sequence  of  pho¬ 
tographs  obtained  by  the  Ranger  spacecraft  as  it  crashed  into  the  lunar  surface. 
This  made  photogram  metric  analysis  dinicult. 


